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Abstract 

Purpose. Radiation therapy is a local treatment aimed at cells in and around a 
tumor. The goal of this study is to develop an algorithmic solution for predicting 
the position of a target in 3D in real time, aiming for the short fixed calibration 
time for each patient at the beginning of the procedure. Accurate predictions of 
lung tumor motion are expected to improve the precision of radiation treatment 
by controlling the position of a couch or a beam in order to compensate for 
respiratory motion during radiation treatment. 

Methods. For developing the algorithmic solution, data mining techniques are 
used. A model form from the family of exponential smoothing is assumed, and 
the model parameters are fitted by minimizing the absolute disposition error, 
and the fluctuations of the prediction signal (jitter). The predictive performance 
is evaluated retrospectively on clinical datasets capturing different behavior (be¬ 
ing quiet, talking, laughing), and validated in real-time on a prototype system 
with respiratory motion imitation. 

Results. An algorithmic solution for respiratory motion prediction (called ExSmi) 
is designed. ExSmi achieves good accuracy of prediction (error 4 — 9 mm/s) 
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with acceptable jitter values (5 — 7 mm/s), as tested on out-of-sample data. The 
datasets, the code for algorithms and the experiments are openly available for 
research purposes on a dedicated website. 

Conclusions. The developed algorithmic solution performs well to be prototyped 
and deployed in applications of radiotherapy. 

Keywords: Respiratory motion compensation, exponential smoothing, 
predictive modeling, real-time 


1. Introduction 

The goal of radiotherapy treatment is to destroy the tumor and at the same 
time prevent the healthy surrounding tissues from being damaged [HElElllllS]. 
Advances in radiotherapy technologies, such as intensity modulated or image 
guided radiotherapy, and stereotactic body radiotherapy, have made highly con¬ 
formal and accurate treatment [6] possible. An important limiting factor to 
the success of tightly conforming dose distributions is the ability to aim the 
radiation beam precisely at the target with minimal positional error. 

Therefore, motion management is one of the most active research and de¬ 
velopment topics in modern radiotherapy, as can be seen from many studies 

n 13 maun]- 

Intrafraction motion (motion of the target during treatment) is usually 
caused by the skeletal muscular, cardiac, gastrointestinal and respiratory sys¬ 
tems, the later being responsible for the most of it. 

The positions of all the organs in the thorax and abdominal regions are af¬ 
fected by respiration of a patient; however, the organs may move in different 
ways and various magnitude. In addition, the tumor itself may be moving along 
with the organs, depending on its location and fixation to the surrounding struc¬ 
tures. The magnitude of the motion highly depends on the location of the tumor 
and also may vary a lot for individual patients. Lung tumors can exhibit up 
to 3 cm motion in the cranio-caudal direction during normal respiration, while 
tumors of other types typically move only a few millimeters or do not move 
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at all m- Movement of lung tumors introduces uncertainty in the position¬ 
ing. To account for this uncertainty the conventional radiation therapy requires 
larger treatment margins, as it is recommended by the International Commis¬ 
sion on Radiation Units and Measurements [HUS]. Extra margins may lead 
to large volumes of healthy tissue being destroyed during radiotherapy treat¬ 
ments. Therefore, while higher doses of radiation therapy may improve survival 
rate, healthy tissue sparing is important to reduce side effects of the organs at 
risk. [5]. 

To cope with this problem various techniques have been considered [2]. Ac¬ 
tive motion compensation umi], such as gated radiotherapy HU [IS], breath- 
hold m or tumor tracking [Miiniiiniiii have been introduced into the clinical 
practice. However, these techniques have limitations, e.g., the total treatment 
time signihcantly increases in case of gated radiotherapy [21) . invasive fiducial 
markers need to be implanted jH] , breath-hold works well only in case of com¬ 
pliant patient. Hence, development of new non-invasive techniques, aimed to 
controlling respiratory motion in radiotherapy, is an important task for the 
modern radiation oncology. Some tracking systems, such as VERO [23], that 
use a beam for positioning and some, like CyberKnife [53] use robotic arm to 
move linac. 

A generic approach to the compensation for respiratory motion is defined 
as follows (following [T]): (i) determine the current position of the tumor from 
an external marker position using a computational technique for relating the 
marker and the tumor [1113Halil]; (ii) predict the next position of the tumor; 
(iii) compensate for the anticipated respiratory motion (e.g. by repositioning the 
beam); and (iv) adapt the dosimetry to the changing configuration of the tumor. 
The current position of a tumor can be determined using external markers [261 
EEaE]. Once the next position of the tumor is predicted, various techniques 
can be used to compensate for the respiratory motion Emm], e.g. shifting the 
patient using a robotic-couch, shifting the beam by repositioning the radiation 
source, redirecting the beam electromagnetically, or changing the aperture of 
the beam. 
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This study focuses on step (ii), i.e., predicting the next position of the tumor 
from the past observations. Prediction is necessary to overcome delays intro¬ 
duced by tracking system latency. For predicting the tumor motion a num¬ 
ber of predictive modeling techniques have been considered [Dull], such as: 
Kalman filters artificial neural networks (ANN) |2[2n], state-based 

probabilistic approaches |20], local regression [19], seasonal autoregressive mod¬ 
els (TVSAR) [3D], autoregressive moving average models (ARMA) |3T], multi- 
step linear methods (MULIN) |32|, and wavelet-based multiscale autoregression 
(wLMS) [33]. 


While most of the existing studies propose new advanced predictive models, 
the complete compensation process itself is understudied. After selecting an ac¬ 
curate predictive modeling technique, it is far from trivial to put it in operation, 
for which a full algorithmic solution is required. Algorithmic solutions should 
include step-by-step instructions for automated data pre-processing, model cal¬ 
ibration for a given patient, adaptation to potential variation in data arrival 
rates, confidence estimation and self-diagnosing mechanisms of the model, and 
potential mode switching (e.g., to a simpler model or no prediction at all). The 
calibration procedure should be done as efficiently as possible in order to mini¬ 
mize preparation time, and maximize utilization of the equipment for treatment. 

This paper proposes a full algorithmic solution for respiratory motion pre¬ 


diction for a selected setup (see sec. 2.11, aiming at minimizing the time for 
model calibration. The predictive performance is evaluated on clinical datasets 
off-line and in real-time on prototype system with respiratory motion imitation. 

Several studies develop controllers for motion compensation [341135j , which 
can be seen as algorithmic solutions, however, their focus is on step (iii), i.e. 
compensating for the anticipated respiratory motion, in Murphy’s classification, 
while our focus is on step (ii), i.e. predicting the next position of the tumor. 

The remainder of the paper is organized as follows. Sect, [^discusses data 
collection (2.11, prediction setting (2.2), performance criteria (|2.3|), prediction 


methods ( |2.4| ), algorithmic solution (2.5) and experimental evaluation 
In Section [^ the performance of the algorithm is evaluated, and in Section [^ 
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Figure 1: General setup for data collection. 

experimental results are discussed. Conclusions and future research directions 
are presented in Section 

2. Materials and Methods 

2.1. Data Collection 

Clinical data is collected using an infrared stereo-camera with 60 Hz internal 
sampling frequency, external markers, HexaPOD evo couch and in-house soft¬ 
ware. Elekta HexaPOD evo RT Systenf] (Elekta AB, Stockholm Sweden) is 
a radiation therapy system setup, with static positioning system iGuide®2.0 
developed by Rubedo Systemt]^ (Rubedo Systems, Kaunas, Lithuania). The 
system was adapted to collect real-time data by Rubedo Systems. 

The radiation treatment system under consideration consists of several com¬ 
ponents: patient setup couch, in this case the HexaPOD coucf0 [36] . radiation 

^http://www.elekta.com/healthcare-professionals/products/elekta-oncology/ 

treatment-techniques/positioning-and-immobilization/hexapod-evo-rt-system.html 
'"http: //rubedo . It/ 

^http://www.elekta.com/assets/Elekta-Oncology/Stereotactic-Radiation-Therapy/ 
case_studies/ 
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Figure 2: Camera position for data collection. 


beam source, usually a medical linear accelerator (linac), tracking device, which 
provides information about the position of the patient and a controller that 
controls the treatment process. Several different control schemes have been 
proposed [371 EH Eg EH 13 . 

Respiratory motion in HexaPOD is measured by an infrared stereo-camera 
(NDI Polaris [20] (NDI (Northern Digital) International, Ontario, Canada)), 
that tracks external markers placed on the body of the patient. We use 1 mm 
spatial resolution, and while the camera can provide up-to 0.25 mm resolution 
under ideal conditions [40] . often it may go up-to 0.6 mm (with 95% confidence 
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interval), therefore 1 mm is the safe choic^ The camera provides position data 
periodically in frames (and frame numbers). The timestamps are computed 
from the sampling frequency of the internal camera, which is 60Hz (a frame in 
each 16.7 ms). Processing delays are negligible (< 1 ms). 

The existing setup (provided by Rubedo systems) is restricted to processing 
every second (2nd) frame, therefore the effective sampling rate is 33.(3), 66.(6) 
or 100 (99.(9)) ms. See general setup and schematic camera position in Figs. 
and[^ (this setup is used for the development and testing of iGuide softwar^. 

As a result, the raw incoming data is not completely equally spaced in time, i.e. 
the time intervals from the second to the sixth or the ninth frame may differ. We 
ensure that the data for prediction is equally spaced by resampling the incoming 
data at a rate that is a multiplier of six frames (which correspond to 100 ms). 
Due to the same reason, the prediction horizon should also be a multiplier of 
six frames. Prediction horizon isselected based on the specifics of setup, where 
we have 100 ms camera communication delay, and we predict future position 
100 ms ahead to compensate velocity of the couch (16 mm/s). 

Ten clinical datasets are used in this study. Each dataset includes 3-dimensional 
observation records with 3 positions per record over time. Each dataset records 
an empty treatment session (no radiation) lasting from 72.617 to 320.05 s. The 
datasets have been collected from 3 healthy males aged 20 to 40. See Table 
for further information about the dataset. 

2.2. Prediction task 

Given is a three-dimensional time series recording the position of an external 
marker over time. The position is given in three coordinates x, y and z in 
millimeters transformed in such a way that min{xi) = Q,min{yi) = 0 and 
min{zi) = 0. Let = {xi,yi, Zi) denote the true position of a marker at time i, 
and let = {x^, , z^) denote the predicted position h steps ahead. When the 

^In case of 4 — 6 mm it could be insufficient. 

®See https://www.youtube.com/watch?v=a4Fqgl6avtA 
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Figure 3: Two predictions giving the same prediction error but different jitters. 

horizon h is clear from the context, it will be omitted from the notation. For 
brevity we index time series by the index of arrival and not by the time-stamp 
of arrival. The index i refers to the number of the current observation in a 
sequence from the start of the reading on the current patient. 

2.3. Performance criteria 

From the operational point of view two performance characteristics are crit¬ 
ical: predictions should be accurate and the predicted signal should fluctuate 
as little as possible (have low jitter [2]). The latter requirement is due to the 
need for the beamer or the couch to move, following the predicted signal, in 
order to compensate respiratory motion. Following sudden rapid movements of 
predictions is impractical and may be infeasible due to mechanical limitations of 
couch or another tracking device. Fig. j^gives two example predictions that lead 
to the same prediction error, but have different jitters. A low jitter is preferred 
from the operational point of view. 

Quantitatively the accuracy of predictions can be measured by a straight line 
distance from the predicted position to the true position in 3-dimensional space 
(3D). For simplicity, distances can be measured in the coordinate units, but 
for interpretability it is better to transform the coordinates and report results 
in standard units of length. This paper reports prediction errors and jitters in 
millimeters. The prediction error at time i is defined as: 

e, = = IId - rill- (1) 

The goal is to minimize the error over a treatment session. Since treatment 
sessions can be of different length, it is practical to look at the mean error over 



a treatment session: 


T 

E = J2e./{TA), ( 2 ) 

i=l 

where T is the duration of the session in number of frames, and A is the time 
interval between two frames. 

The jumpiness or jitter |14j can be measured as the distance the prediction 
signal travels per time step: 

j^ = \/{Xi - + [y, - + {Zi - 

= (3) 

For the units (mm) to be interpretable and comparable to the error, in the 
experiments we will report average jitter and average error per second (A = 0.1). 

The goal is to minimize jitter over a treatment session. Since treatment 
sessions can be of different length, it is practical to look at the mean jitter over 
a treatment session: 

T 

^ = EM(r-l)A), (4) 

i=2 

where T is the duration of the session in number of frames, and A is the time 
interval between two frames. 

Note that jitter is minimized when fi = fi-i for all i G [2,T], i.e. the 
prediction is constant. However, in this case no compensation for respiratory 
movement is possible. In practice, a system aims at compensating for respiratory 
movement, it needs to hnd a balance between error and jitter. 

2.4- Predictive modeling techniques 

We are aiming at developing an algorithmic procedure for real time predic¬ 
tion of respiratory motion. Such an algorithm takes a base model as input, and 
determines when the model should be calibrated, when the actual operation can 
start, and how to switch between alternative models of different complexity. 

For predicting the tumor motion a range of predictive models have been 
considered mmia, as discussed in the introduction. Our main qualitative 
criteria for choosing an existing technique for the algorithm are as follows. 
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1. Models need to be fast to calibrate (up to 30 — 60 sec) for every next 
patient, since waiting time is costly. The number of model design and 
calibration parameters should be minimal. 

2. Models should be able to adapt to changes in respiration rhythm and drifts 
of the tumor during a session. 

3. Models and prediction decisions should be transparent (how predictions 
are made), so that the technique can be trusted by medical specialists. 

4. Models which are simple to implement on any treatment hardware with 
minimal usage of external tools are preferable to minimize risks of software 
errors and dependencies. 

Table provides a summary of considered base models and our assessment 
against the four qualitative criteria. The main limitation of state probabilistic 
methods (such as Kalman filters or Hidden Markov models) and autoregressive 
models (such as autoregressive moving average models, regression models fit¬ 
ted using least squares procedure) is that they require relatively large training 
sample for model calibration before it can be used for predictions, and we are 
looking for very fast and robust models. More advanced machine learning mod¬ 
els (such as neural networks or support vector machines) require even larger 
training sample sizes, and in addition, the resulting models are so called ’’black 
box” models, where it is extremely difficult to trace how the predictions are 
made. Therefore, given the focus of our study on fast, interpretable, adaptive 
and transparent prediction making, we resort to extrapolation and exponential 
smoothing techniques for our algorithmic solution. The next subsections discuss 
these two types of techniques in detail. 

2.4-1- Extrapolation methods 

Extrapolation methods, predict based on the most recent observations. They 
do not require any calibration and minimum or none parameter settings. These 
methods have very short memory of the past data and hence are inherently 
adaptive to changes in respiration rhythm or tumor drifts, they are very trans¬ 
parent (easy to explain to a non-specialist) and very simple to implement. 
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Persistent prediction. (PP) is the simplest predictor, which predicts that the 
next signal will be the same as the last observed. No parameters are required. 

xt+h = xt, (5) 


here t is the time index and h is the prediction horizon. 

Persistent prediction can be considered as a baseline for compensation for 
respiratory motion. It does not predict pro-actively, but only follows past ob¬ 
servations. 

Linear extrapolation. (LE) assumes that the signal will maintain the same ve¬ 
locity and direction as last observed. No parameters are required. 

Xt+h = Xt + {xt-Xt-h)- ( 6 ) 

Multi-step linear prediction (MULIN) |32j is a generalization over linear ex¬ 
trapolation, it takes into account acceleration of the signal of different order. 
Since the extrapolations may become unstable if the signal is noisy, MULIN uses 
exponential smoothing moving average of the predictions instead of outputting 
only the latest prediction. 

Xt+h = a^Xt-\-'Y2Hxt,h)'"'^ -\-{l-a)xt+h-i (7) 

5{xt,h)^ = xt-xt-h ( 8 ) 

5{xt,hy^^ = 6{xt,hy - S{xt-h,hy (9) 

where k G [1,2,...] and a G (0,1) are user specified parameters. In this paper 
we experiment with the second order MULIN. We used the default parameter 
settings supplied in the implementation made available by the author^ 

Exponential smoothing 

Exponential smoothing is a type of moving average, where the importance 
of the past observations decreases exponentially. Exponential smoothing is not 


^http://www.rob.uni-luebeck.de/~ernst/dateien/mulin/mulin.m 
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parameter intensive, the only parameter is the speed with which old observations 
are forgotten. Exponential smoothing does not require model calibration for 
every new patient, it can predict immediately after the start, but a short warm¬ 
up period is advisable. Just like extrapolation methods, exponential smoothing 
is inherently adaptive, similarly transparent, and straightforward to implement. 

Simple exponential smoothing. (ESI) makes prediction as the exponentially 
weighted moving average of the previous observations. 

xt+h = axt + {1 - a)xt-i, ( 10 ) 

for any horizon h. Here a G (0,1) is a user defined parameter. If the forgetting 
factor a is small, then forecasting will have a long memory. If a is close to 
one, then forgetting will be fast, a = 1 would mean that we predict the next 
observation to be the same as the last (PP). a = 0 would give a constant 
prediction (zero jitter). ESI is equivalent to autoregressive integrated moving 
average model [H] ARIMA(0,I,I). 

Simple exponential smoothing does not do well when there is a trend in the 
data. 

Double exponential smoothing. (ES2) takes trends into account. 


Xt+h = h + hbt 

(11) 

h = Q;a;t-|-(1 — Q;)(^t-i + ^t-i) 

(12) 

t't = P{lt — h-i) + (1 ~ P)bt-i 

(13) 


Here a G (0,1) and /3 G (0,1) are user specified parameters. Initialization: 
lo = Xq, bo = 0. ES2 is equivalent to ARIMA(0,2,2). 

In case of double exponential smoothing for respiratory motion prediction 
breath cycle will be modeled as short term trends. 

The main limitation of this approach is that the prediction will systemati¬ 
cally overshoot when the direction of the signal reverses. 
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Holt-Winters exponential smoothing. , or triple exponential smoothing (ESS), 
is often used for short term forecasting of seasonal time series [42], as it can 
handle trends and seasonality. Seasonality means that the signal is periodic 
with a period p. We consider ESS model with additive seasonality component 
(based on US]). 


^t+h 

= It hbt St-p+h 

(14) 

It 

= oi{xt — St-p) -b (1 — a){lt-i -b bt-i) 

(15) 

bt 

= /3(4 ~ ^t-i) + (1 ~/3)&t-i 

(16) 

St 

= 7{xt - It-i - bt-i)(1 - 7)st-p 

(17) 


Here a S (0,1), ,0 S (0,1) and 7 G (0,1) are user specified parameters. Initial¬ 
ization: Iq = Xo, 60 = 0, So, ... , St-p = 1. 

The original ESS requires the period to be known and fixed during the 
model operation. The period of a respiratory signal, however, varies even for 
a single patient, as respiration may become more frequent or slow down over 
time, the air intake may be delayed due to talking or coughing. We make a 
stabilizing modification to ESS, we use the initial level in estimation of the 
seasonal component instead of moving average of the level: 

St = 7(xt - lo) + (1 - 7)st-p- (18) 

We suggest using the parameter values listed in Table which have been 
found during initial experiments on the training parts of a couple of traces. 
The testing part of the traces on which the accuracies are reported, was never 
used for estimating these parameters. To minimize the chance of overfitting the 
training data the parameters are fixed for all the traces. 

We suggest using a fast forgetting for the level (having in mind potential 
bias of the model and potential drifts), keeping it within a recommended [43] 
restriction 0 < a -I- 7 < 1. The role of the trend component is to estimate 
long term changes in the average signal level, thus the memory should be long, 
thus for ESS /3 should be low. Since ES2 has no seasonal component, the trend 
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component plays the role of seasonal adjustment, thus /3 needs to be higher than 
in ESS, but not too high, since in such a case overshooting at turning points may 
be too large. Since we know that ES2 is biased (data contains seasonality, but 
we approximate it by the trend component), we need to have a fast forgetting 
not to propagate model bias therefore a should be high. 

£.5. Prediction procedure ExSmi 

We propose the following procedure for predicting respiratory motion, called 
ExSmi, summarized in Algorithmic ExSmi includes online preprocessing out¬ 
lier removajC (condition on line 9), online model calibration and switch pre¬ 
diction phase (line 11), a switching mechanism between the main model and a 
simple, but more robust baseline (line 18), which is based on the most recent 
performance, taking into account two quantitative criteria - prediction error and 
jitter (line 18). Linear extrapolation method is considered as a baseline S, and 
exponential smoothing is considered as the main predictive model (L). 

At the time of model switch (line 18) both models are well warmed up, and 
estimates of the most recent performance are available. We select the model 
demonstrating the lowest recent prediction error and jitter of the two and apply 
a fading factor a to the running estimates of the performance to ensure that the 
most recent performance is accounted with more weight, while considerng not- 
so-recent performance history with lower weight. IT helps to minimize the risk 
of sudden jumps in prediction error or jitter, when the predictors are switched. 
In the next section we experimentally analyze the performance of the proposed 
approach. 

A recent study [H] aims at classifying the patients into predictable and 
unpredictable, in order to decide whether motion compensation should be used 
at all, which is conceptually similar to our approach, but there are several key 
differences. While the authors consider whether motion compensation should 

^The threshold has been selected based on speed of the coach movement. It is not possible 
that the couch moves that fast as to produce 1 cm difference between points. 
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be used or not, we do not question the applicability of motion compensation, 
but dynamically switch between models of different complexity depending on 
the noisiness of the signal. Moreover, their approach is to decide regarding each 
patient before commencing the treatment, while our approach is intended for 
real time, and the outcome during each new treatment session may be different. 

2.6. Protocol for experimental analysis 

We experimentally analyze the performance of the base models selected in 
Section [2^ and the proposed algorithmic solution in the following settings. New 
observations arrive every 100 ms and the required prediction horizon is 200 ms 
ahead (h = 2). The warm up period is 30 sec, which is 300 samples {w = 300). 
Prediction errors and jitters are reported as averages from observation 301 until 
the end of the treatment session. We first test the prediction methods stand 
alone, and then test a selected prediction method inside the proposed algorithm. 

All the experiments are performed using in-house produced MATLAB@ 
code, available at http://datasets.bpti.lt/radiotherapy 

3. Results 

Our experimental analysis consists of two parts: firstly, we experimentally 
evaluate the performance of alternative prediction methods in terms of predic¬ 
tion error and jitter, and secondly, we experimentally analyze the performance 
of the proposed algorithmic solution. 

3.1. Predictive performance of base models 

Figure [^depicts prediction errors and jitters of the base models, selected in 
Section [ 2 )^ On the left plot each dot is one time series (recall that each dataset 
includes three positions, that is, three time series). We can see from the left plot 
that the selected models provide a variety of errors and jitters, indicating that 
some of the signals are more difficult to predict than the others. However, dots of 
the same color (the same base model) appear in elongated clusters, suggesting 
that there may be a trade-off between error and jitter achieved by different 
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Figure 4: Predictive performance of alternative base models, (left) Each dot is one time series, 
(right) average performance. 

models, that is, a gain in error increases jitter and the other way around. Some 
models may be better at minimizing error, others at keeping the jitter low. 
The right plot presents the average overall time series for each model. We see 
that ESI, PP, ESS and LE demonstrate nearly a linear trade-off between jitter 
and error with ESI showing the lowest jitter and LE showing the lowest error. 
MULIN demonstrates a reasonable error, but the jitter is much worse than that 
of the other models. ES2 achieves nearly the same error as LE, but has a lower 
jitter, therefore, we select ES2 as the primary base model for our algorithmic 
solution. The performance of a constant prediction, which achieves zero jitter, 
is not plotted since the error (51mm/s) is too far off the scale of the plot. 

3.2. Visual analysis 

Eigurej^ plots the predictions of the compared methods on a snapshot of the 
first coordinate from the experiment 201205181211-UAC-1-N-320-6. We can see 
that PP and ESI have a regular delay in predictions, LE and ES2 overshoot at 
peaks, MULIN and ESS follow the signal reasonably well, but MULIN is too 
spiky (high jitter) and ESS occasionally makes sudden errors. Based on this 
visual analysis ES2 or LE are preferred methods. Figure compares the 
jitters produced by different methods. Each line shows how much the beam 
would need to travel in 10 seconds if the predictions were followed. We see that 
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Figure 5: Performance on a snapshot of experiment 201205181211-UAC-1-N-320-6. 
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Figure 6: How much the couch would need to travel in 10 s if the predicted signals are followed 
(including both jitter, and respiratory motion), averaged over experiment 201205181211-UAC- 
l-N-320-6. 


all methods are comparable in terms of jitter except for MULIN, which produces 
substantially larger jitter. 

Next we look at the scatter of predictions in space from a patient’s perspec¬ 
tive. Assuming that the bed can perfectly track the predictions, we plot where 
the beam will hit in 2D with respect to the true target. For that we subtract 
the true signal from the prediction, this way the true target is always (0,0). 
Intuitively, in order not to cause unnecessary harm to the patient, deviations 
of predictions from the target (0, 0) should be as small as possible and there 
should be no far outliers. Moreover, the errors should be distributed around the 
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Figure 7: Scatters of predictions on the testing range of 201205181211-UAC-1-N-320-6. 


target (0, 0) as evenly as possible, not concentrated in one or a few spots. 

Fig. [Tjplots the scatters of predictions for the same experiment 201205181211- 
UAC-l-N-320-6 ) (resampled 201205181211-UAC-1-N-320) in 2D (x and z co¬ 
ordinates). We see that all the six methods produce predictions that are rea¬ 
sonably close to the true target, as compared to no compensation. However, 
PP and LEI have the strongest tendency to make concentrated errors, meaning 
that particular two spots on the upper right and lower left sides from the target 
may be burned due to prediction latency. We would like to notice that in this 
research tumor is treated as a point (centroid) representing a 3D volume. 

3.3. Performance of the proposed algorithm 

We investigate the performance of ExSmi algorithm with the second order 
exponential smoothing ExSmi(ES 2), which showed the most promising perfor¬ 
mance in the previous experiment. We compare the performance of the al¬ 
gorithm ExSmi(ES 2) with applying ES2 and a naive persistent prediction PP 
stand alone. Figure plots the errors and jitters on all experiments, one dot 
represents one dataset and the numerical results are provided in an on-line ap¬ 
pendix]^ We can see from the plot that ExSmi(ES 2) has advantage over simple 


® On-line appendix http: //datasets .bpti. lt/wp-content/uploads/2013/12/ 

ExSmi-OnlineAppendix.pdf 
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All datasets 



error (mm/s) 


• PP 

• ES2 

• ExSmi(ES2) 


Figure 8: Prediction accuracy. 


ES2 in situations where overall error and jitter are quite high, i.e. in extremely 
unpredictable cases. This performance supports the intuition, that when an 
intelligent method cannot do well, it makes sense switching to a robust baseline 
predictor. 

Next we analyze this phenomenon in more detail. Table provides average 
errors and jitters for the experiments divided into two groups: (1) difficult to 
predict identified by high prediction error (> 8 mm/s) and easy to predict 
identified by lower prediction error (< 8 mm/s). We see that, indeed, for the 
difficult to predict cases the algorithm provides a better balance between error 
and jitter, while it does not disturb much the easier to predict cases. 

Our dataset includes signals with different activities (such as laughing or 
talking). Next we analyze the performance of ExSmi at different activities. 
Table presents the results. We can see that normal position demonstrates 
the lowest overall error and jitter, as it could be expected, since the patient 
stays still. Prediction in laughing and talking conditions is, hence, more diffi¬ 
cult. The proposed ExSmi performs nearly the same as the base model ES2 
in normal/other conditions; however, ExSmi consistently performs the best in 
other than normal conditions, which is a desired feature of our solution. We 
have implemented outlier control and predictor switch mechanisms so that the 
predictions stay robust in difficult situations, and these experimental results 
support that. 
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full treatment cycle, AP directions (Z coord) 


snapshot 


snapshot 



Figure 9: Performance on a snapshot of experiment 201205101536-LAC-l-LT-142-6. 


Finally, Figure|^plots an example of predictions by the algorithm ExSmi(ES 2) 
and ES2 stand alone on a difficult to predict case. We see that when the true 
signal suddenly starts to jump ES2 largely overshoots. This is because ES2 
takes into account the velocity of the signal, observing one sudden jump in the 
signal level leading to extrapolation of this pattern, i.e. predicting that the 
signal will jump further with a similar speed. In such a case when the signal 
is noisy a naive persistent predictor proves to be more accurate. The proposed 
algorithm combines ES2 and PP and takes advantage of both. 

4. Discussion 

ExSmi, PP and LE approaches have been implemented in a prototype 
Rubedo system including a HexaPOD couch, and an infrared stereocamera (NDI 
Polaris). This prototype validation has confirmed our experimental results, and 
several additional observations regarding the performance have been made. 

1. HexaPOD couch is quite sensitive to larger speed and direction changes 
and jitter, i.e. the device starts vibrating. Currently, the problem is 
solved by putting an independent restriction on velocity changes. As a 
future work, it would be interesting to consider such constraints as part 
of the prediction algorithm. 

2. ExSmi(ES 3) implementation seems to be over-sensitive to periodicity changes, 
while the period of respiratory motion typically is changing all the time. 
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That explains why the best results have been achieved by ExSmi(ES 2) 2 
and LE with anti-vibration means. 

In this paper only prediction of tumor motion and its compensation is in¬ 
vestigated, but in case of external markers motion of tumor should be predicted 
from a motion of an external marker, e.g. |451146j . which would induce addi¬ 
tional error. This technique requires fixing markers near tumor. However they 
should be out of the beams’ way. The couch used in the setup is not constructed 
to compensate motion, therefore it is not clear how long and how well it would 
operate over the extended period of time in such a mode. 

Experiments were performed on motion recorded using external markers 
under an assumption that tumors move in a similar fashion. Therefore, further 
investigation with tumor motion could be useful. 

The important question, which we did not answer in this paper, is how much 
would a prediction would correct a clinical misalignment of the target? It could 
be, that linac, MLC, immobilization devices and, especially, a live patient are 
contribute more the the overall error, while precisions of the most of the existing 
predictors is sufficient. It is out of scope of this paper, but such analysis could 
be very interesting. 

In summary, the prototype implementation has demonstrated a promising 
performance, confirmed our experimental findings, and indicated an interesting 
direction for future investigation. 

5. Conclusions 

The study investigated prediction models and developed an algorithmic so¬ 
lution ExSmi for predicting the position of a target in 3D in real time. ExSmi 
demonstrated good performance, measured by the prediction accuracy and the 
jitter of the prediction signal. The developed algorithmic solution performs well 
to be prototyped deployed in radiation therapy applications. 

This study has opened several interesting and important directions for fu¬ 
ture research. The first direction is to extend the algorithmic solution ExSmi 
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to take into account technical characteristics of the equipment, for instance, the 
maximum possible velocity change of the treatment couch. While this study 
treated each respiratory signal as an independent observation, the second in¬ 
teresting and important direction for extension would be to consider multiple 
signals from different locations simultaneously. Taking into account such char¬ 
acteristics it is expected to further improve the precision of treatment. 
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Appendix A. Jitter per time spacing and relative error 

Table [576] presents jitter and error of Figure [^relative to PP. 

Table [X?7| presents jitter and error of Table relative to PP. 
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Table 1: Experimental Data, where SI (superior-inferior) positions are (L-lower, U-upper), 
body positions (A-abdomen, C-chest and LR (L-left, C-center, R-right), and different states 
(T-talking, N-normal, O-other (other type of motion), L-laughing); directions: SI (superior- 
inferior), LR (left-right), AP (anterior-posterior) 


Signals 

Max Range 

SI, LR, AP (mm) 

Duration 

(s) 

Frames 

Experimental setting 

201205101519- 

LACUACUCC-3- 

T-222 

19, 

4, 23 

222.00 

6658 

lower abdomen cen¬ 
ter, upper abdomen 

center, upper chest 

center, talking 

201205101522- 

LACUACUCC-3- 

N-138 

6, 

3, 20 

138.00 

4148 

lower abdomen cen¬ 
ter, upper abdomen 

center, upper chest 

center, normal 

201205101534- 

LACUACUCC-3- 

NO-130 

9, 

4, 20 

130.00 

3883 

lower abdomen 

center, upper 

abdomen center, 

upper chest center, 

normal, other 

201205101536- 

LACUACUCC-3- 

LT-142 

29, 

14, 31 

142.00 

4267 

lower abdomen 

center, upper 

abdomen center, 

upper chest center, 

laughing, talking 

201205101541- 

LACUACUCC-3- 

N-130 

6, 

2, 17 

131.00 

3919 

lower abdomen cen¬ 
ter, upper abdomen 

center, upper chest 

center, normal 

201205111055- 

LACLARUAR-3- 

N-117 

6, 

4, 18 

117.00 

3513 

lower abdomen cen¬ 
ter, lower abdomen 

right, upper ab¬ 
domen right, nor¬ 
mal 

201205111057- 

LACLARUAR-3- 

0-72 

40, 

10, 45 

72.62 

2178 

lower abdomen cen¬ 
ter, lower abdomen 

right, upper ab¬ 
domen right, other 

201205181211- 

LACUACUCC-3- 

N-320 

12, 

4, 31 

320.05 

9593 

lower abdomen cen¬ 
ter, upper abdomen 

center, upper chest 

center, normal 

201205181220- 

LACUACUCC-3- 

N-306 

20, 

5, 36 

306.00 

9176 

lower abdomen cen¬ 
ter, upper abdomen 

center, upper chest 

center, normal 
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Table 2: Qualitative assessment and selection of base models. 


Technique 

Fast to 

calibrate 

Adaptive 

Transparent 

Simple to 

implement 

Select 

Extrapolation methods 

yes 

yes 

yes 

yes 


Exponential smoothing 

yes 

yes 

yes 

yes 


State-based probabilistic 

no 

yes/no 

yes 

yes 


Autoregressive models 

no 

yes/no 

yes 

yes/no 


Neural networks 

no 

yes/no 

no 

no 


Support vector machines 

no 

no 

no 

no 



Table 3: Recommended parameters for exponential smoothing. 


Model 

level 

a 

trend 

/3 

seasonal 

7 

respiratory rate 

P 

ESI 

0.7 




ES2 

0.7 

0.6 



ESS 

0.7 

0.3 

0.3 

5.5 sec 


Table 4: Average performance on difficult and easy to predict cases. 


Group 

Measure (mm) 

PP 

ES2 

ExSmi)!?^^) 

Difficult 

average error 

10.8 

9.6 

9.1 

E 

(std.) 

(2.1) 

(1.5) 

(1.3) 

> 

average jitter 

6.0 

8.8 

7.0 

8 mm/s 

(std.) 

(1.1) 

(1.5) 

(1.2) 


error + jitter 

16.8 

18.4 

16.1 

Easy 

average error 

7.5 

4.1 

4.2 

E 

(std.) 

(3.7) 

(1.6) 

(1.6) 

< 

average jitter 

4.1 

5.0 

4.8 

8 mm/s 

(std.) 

(1.8) 

(1.9) 

(2.0) 


error + jitter 

11.5 

9.0 

9.0 
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ALGORITHM 1: Predict respiratory signal h steps ahead 


1: incoming observations r = (x,y,z) 

2: predictive model form L with design parameters 6 
3: prediction horizon h, warm-up w (recommended ~ 30 s) 
4: decay for measuring recent error d G (0,1) (rec. d — 0.1) 


5: Initialize model Lq (See Sec. 2.4 for recommenations) 

6: Initialize error and jitter counts Eq — 0, Eq = 0, = 0, = 0 

7: for t - 1 — 2,..., / do /*from the start to the end of treatment*/ 

8: receive the latest observation rt 

9: if \\rt — rt_i|| < 1 cm then 

10: update model Lt = f{Lt-i, {x, y, z)t) 

11: if t < w then /*if warmup is over make predictions*/ 

12: make prediction with Lt'. 

13: make baseline prediction = n + [rt — rt-h) 

14: error Ef = d* error{rt,rt) + (1 — d)Ef_i [Eq. Q] 

15: error Ef = d* error{rt-h, rt) + (1 — d)Ef_i [Eq. Q] 

16: jitter j/ = d* jitter{rt, ft-i) + (1 — d)j/_i [Eq. (j^[ 

17: jitter = d * jitter{rt-h,rt-h-i) + (1 - d)j(Li [Eq. ([^] 

18: if {Ef + j/) > {E/ + jf) then /*L performs well*/ 

19: final prediction by the main model ft+h = ft+h 

20: else 

21: final prediction by baselineft+h = 

22: end if 

23: end if 

24: else/*rt is an outlier, ignore*/ 

25: if t < w then 


26: predict -ft+h = ft+h-i 

27: set Et = E/_t, E/ = E(L„ J/ = j/_,, jf = 

28: end if 

29: end if 

30: adjnst the beamer /*out of the scope of this paper*/ 

31: end for 
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Table 5: Average performance with different acitivities. 


Group 

Measure (mm/s) 

PP 

ES2 

ExSMi(i?5.2) 

Normal 

average error 

7.0 

3.9 

4.0 


average jitter 

3.9 

4.7 

4.6 


error + jitter 

10.9 

8.6 

8.6 

Normal/ 

average error 

8.1 

4.0 

4.0 

other 

average jitter 

4.2 

5.0 

5.0 


error + jitter 

12.3 

9.0 

9.0 

Talking 

average error 

10.2 

4.6 

4.8 


average jitter 

5.3 

6.2 

6.1 


error + jitter 

15.6 

10.8 

10.9 

Talking 

average error 

7.6 

6.1 

6.1 


average jitter 

4.2 

5.9 

5.1 


error + jitter 

11.8 

12.0 

11.2 

Laughing/ 

average error 

10.8 

10.1 

9.4 

talking 

average jitter 

6.0 

9.1 

7.0 


error + jitter 

16.8 

19.2 

16.4 


Table A.6: Predictive performance of the base models relative to persistent prediction (PP). 



PP 

LE 

MULIN 

ESI 

ES2 

ESS 

relative error 

1.00 

0.60 

0.74 

1.17 

0.62 

0.78 

relative jitter 

1.00 

1.45 

2.06 

0.92 

1.27 

1.22 
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Table A.7: Relative performance on difficult to predict and easy to predict cases. 



Measure 

PP 

ES2 

ExSmi(ES2) 

Difficult 

relative error 

1.00 

0.89 

0.84 

E>8 

relative jitter 

1.00 

1.47 

1.17 

mm/s 

error + jitter 

1.00 

1.09 

0.96 

Easy 

relative error 

1.00 

0.55 

0.56 

E<8 

relative jitter 

1.00 

1.21 

1.18 

mm/s 

error + jitter 

1.00 

0.78 

0.78 
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